/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      fu;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 0 0 0 0];

internalField   #codeStream
{
    codeInclude
    #{
        #include "fvCFD.H"
        #include "constants.H"
    #};

    codeOptions
    #{
        -I$(LIB_SRC)/OpenFOAM/lnInclude \
        -I$(LIB_SRC)/finiteVolume/lnInclude \
        -I$(LIB_SRC)/meshTools/lnInclude
    #};

    code
    #{
        const IOdictionary& d = static_cast<const IOdictionary&>(dict);
        const fvMesh& mesh = refCast<const fvMesh>(d.db());

        const vectorField& points = mesh.C();
        scalarField f(mesh.nCells());

        scalar H2molf = 0.25;

        scalar MH2 = 1.00784*2.0;
        scalar MN2 = 14.0067*2.0;
        scalar MO2 = 15.9988*2.0;
        forAll (points, celli)
        {
            scalar y = points[celli].y() - 0.03;
            scalar _10H2 =
              - 3.50957e+04*pow(y,5)
              - 3.26065e+04*pow(y,4)
              - 5.25765e+03*pow(y,3)
              - 2.62149e+02*pow(y,2)
              - 3.84660e-01*y
              + 1.99094e-01;
            scalar _20H2 =
                1.84717e+05*pow(y,5)
              - 1.40719e+04*pow(y,4)
              - 6.78365e+03*pow(y,3)
              - 4.17966e+02*pow(y,2)
              - 6.09813e-01*y
              + 3.77362e-01;
            scalar _30H2 =
                3.70457e+05*pow(y,5)
              + 1.89256e+04*pow(y,4)
              - 5.74215e+03*pow(y,3)
              - 4.60436e+02*pow(y,2)
              - 4.98623e-01*y
              + 5.26445e-01;
            scalar _40H2 =
                5.69961e+05*pow(y,5)
              + 6.66643e+04*pow(y,4)
              - 2.62118e+03*pow(y,3)
              - 4.11878e+02*pow(y,2)
              - 1.79756e-01*y
              + 6.48390e-01;
            scalar xH2;
            if(H2molf < 0.10)
            {
                xH2 = H2molf/0.10*_10H2;
            }
            else if (H2molf < 0.20)
            {
                xH2 = _10H2 + (H2molf - 0.10)/0.10*(_20H2 - _10H2);
            }
            else if (H2molf < 0.30)
            {
                xH2 = _20H2 + (H2molf - 0.20)/0.10*(_30H2-_20H2);
            }
            else if (H2molf <= 0.40)
            {
                xH2 = _30H2 + (H2molf - 0.30)/0.10*(_40H2 - _30H2);
            }
            else
            {
                Info<<"\n Error: molar H2 fraction above 0.40 \n";
            }

//             scalar xH2 =
//                 coeffs[0]
//               + y*coeffs[1]
//               + y*y*coeffs[2]
//               + y*y*y*coeffs[3]
//               + y*y*y*y*coeffs[4]
//               + y*y*y*y*y*coeffs[5];
            scalar xO2 = (1.0 - xH2)*0.21;
            scalar xN2 = 1.0 - xH2 - xO2;

            f[celli] = xH2*MH2/(xN2*MN2 + xH2*MH2 + xO2*MO2);
        }
        os  << "nonuniform " << f;
    #};
};

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    walls
    {
        type            zeroGradient;
    }
    defaultFaces
    {
        type            empty;
    }
}


// ************************************************************************* //
